Решение уравнения Фолкнера-Скэн в среде Mathcad

Глава из книги «Решение дифференциальных уравнений в Mathcad»,

готовящейся к публикации в издательстве МЭИ

The solution of the Falkner-Scan equation with Mathcad

Alexander P. Solodov

Table of contents:

Introduction

The mathematical formulation

Reduction of a boundary-value problem to an initial problem by the  sbval method

Solution of an initial problem by the rkfixed method

Imaging of a flow field

Boundary layer on a permeable surface.

Thermal boundary layer equation

Troubles at the solution of a boundary-value problem by means of Odesolve function

Conclusion

Literature

 

The nonlinear differential equation of the third order known as the Falkner-Scan equation

at the relevant boundary conditions

describes the class of so-called similar laminar boundary layer flows  on a permeable surface (see figure) and at varied velocity (varied pressure) in an exterior stream.

Dependent variable f  is a dimensionless stream function, independent variable h - dimensionless distance from a wall. A longitudinal component of velocity and tangential stress are defined as first and second derivative from f. The solution of a boundary-value problem serves the base to the analysis such important problems, as   

q              a curvilinear contours flow (blades of gas turbines or wings of the plane), 

q              intensive condensation or evaporation on interphase,

q              a mass transfer in a boundary layer at intensive catalytic reaction on a wall,

q              processes with the specially organized blowing or suction through a permeable wall.

Numerical parameter fw  in boundary conditions on a wall sets the mass cross-flow on a surface. The positive values determine flows with suction, negative - with blowing. The zero value corresponds to the flow on a impermeable surface. Numerical parameter b (positive or negative) in the Falkner-Scan equation sets a degree of acceleration or deceleration of an exterior stream. Below mainly the flow with zero value of this parameter  will be surveyed.

The Falkner-Scan equation, together with the thermal boundary layer equation, similar on structure, forms the theoretical base of a convective heat mass transfer. Both mentioned ordinary differential equations are received from more blanket partial differential equations of impulse and energy transport, on the basis of the following two fundamental ideas:  asymptotic conversion to the boundary layer equations at the high Reynolds numbers and  then similarity transformation.

However even being result of essential simplifications, the Falkner-Scan equation remains by composite mathematical object, owing to nonlinearity and expedient of the boundary conditions assignment. We shall show further, that the mathematical package Mathcad  gives effective tools for the solution of such uneasy problems.

Main objective will be achievement of understanding, how the viscous stream  and temperature  field  near to wall  is shaped  and how it is possible,  to affect  the flow structure  by means of exterior actions, such as a blowing or suction. We shall try to achieve this purpose in such a manner that we shall change driving parameters and to make the results visual by means of the distribution diagrams and of vector fields pictures (see figure).

Parallel with it, the reader can become proficient in Mathcad technique, namely,  in  numerical integration of two-point boundary problems, in plotting, in development of the user function with  the built-in  programming language etc.

 

 

 

 

 

 

Решение уравнения Фолкнера-Скэн в среде Mathcad

А.П. Солодов

Введение

Математическая формулировка

Сведение краевой задачи к начальной задаче методом   sbval

Решение начальной задачи методом rkfixed

Построение поля течения

Пограничный слой на проницаемой поверхности.

Уравнение теплового пограничного слоя

Неприятности при решении краевой задачи посредством функции Odesolve

Заключение

Литература

Введение

Нелинейное дифференциальное уравнение третьего порядка, известное как уравнение Фолкнера-Скэн [1],

( 1)

при граничных условиях

( 2)

описывает класс так называемых автомодельных (подобных)  ламинарных течений в пограничном слое на проницаемой поверхности (Рис. 1) и при изменяющейся скорости (изменяющемся давлении) во внешнем потоке. Зависимая переменная f есть безразмерная функция тока, независимая переменная h - безразмерное расстояние от стенки. Продольная составляющая скорости и касательное напряжение трение определяются как первая и вторая производные от f. Исследование решений краевой задачи служит фундаментом  для анализа таких важных для практики задач, как  

 

 

Рис. 1 Схема обтекания проницаемой поверхности.

 

Числовой параметр fw  в граничном условии на стенке задает величину поперечного потока массы на поверхности. Положительные значения определяют течения с отсосом, отрицательные – со вдувом (Рис. 1) через проницаемую поверхность. Нулевое значение соответствует течению на непроницаемой поверхности.

Числовой параметр b (положительный или отрицательный)  в уравнении Фолкнера-Скэн задает степень ускорения или замедления внешнего потока. Ниже будут в основном рассмотрены течения с нулевым значением этого параметра, т.е. течения, не подверженные воздействию сильного продольного градиента давления.

Уравнение Фолкнера-Скэн ( 1), совместно с аналогичным по структуре уравнением  ( 17) для теплового пограничного слоя, образует  теоретический фундамент конвективного тепломассообмена. Оба упомянутых обыкновенных дифференциальных уравнения получены из более общих уравнений переноса импульса и энергии в частных производных на основе следующих двух фундаментальных идей.

Первая из них использует то обстоятельство, что   большинство реальных течений происходит при больших числах Рейнольдса, благодаря чему   становится возможным асимптотический переход  к радикально более простым уравнениям пограничного слоя, которые, однако, остаются еще уравнениями в частных производных.

Вторая идея позволяет перейти от уравнений в частных производных к обыкновенным дифференциальным уравнениям благодаря преобразованию подобия - такой замене переменных, что две независимые переменные x и у комбинируются в единственную независимую переменную h=y/ d(x), где d(x) – толщина пограничного слоя (Рис. 1). Другими словами, переменная подобия h есть расстояние от стенки, измеренное в масштабе толщины пограничного слоя.

Однако даже будучи результатом существенных упрощений, уравнение Фолкнера-Скэн остается сложным математическим объектом, вследствие нелинейности и способа задания краевых условий [2,3].  Мы покажем далее, что математический пакет Mathcad предоставляет эффективные средства для решения таких непростых задач.

Главной целью будет достижение понимания того, как формируется поток вязкой жидкости и поле температур вблизи обтекаемой поверхности и как можно повлиять на структуру течения посредством внешних воздействий, таких как вдув или отсос. Мы постараемся достичь этой цели, варьируя в расчетах управляющие параметры и делая наглядными, видимыми результаты -  с помощью графиков распределений, картин векторных полей и линий тока. Параллельно читатель может овладеть техникой работы в среде Mathcad, а именно, численным интегрированием системы обыкновенных дифференциальных уравнений с краевыми (двухточечными) условиями, построением диаграмм различного типа, разработкой пользовательской функции с помощью встроенного языка программирования и т.д.

Математическая формулировка

Поскольку основной интерес представляют профиль скорости и трение, т.е. первая и вторая производные от решения f, целесообразно представить дифференциальное уравнение  третьего порядка ( 1) в виде системы трех дифференциальных уравнений первого порядка. Искомую функцию, ее первую и вторую производную будем рассматривать как компоненты вектор-функции  (F0, F1, F2) . Проделав следующие очевидные по смыслу замены

получим систему уравнений в векторной записи:

( 3)

где зависимая переменная и правая часть определены как вектор-функции:

 

( 4)

Поскольку  полученная система дифференциальных уравнений нелинейна (содержит квадратичные члены F0×F2 и F12), потребуется численная процедура интегрирования.

Краевые условия ( 2)   следует переписать  в новых обозначениях следующим образом:

 ( 5)

Первое из этих условий указывает значение безразмерной   функции тока  на стенке (при h=0). Это варьируемый параметр, посредством которого задается интенсивность вдува или отсоса.

Второе условие задает нулевое значение продольной составляющей скорости на стенке (при h=0), что соответствует фундаментальному условию прилипания.

Третье условие определяет продольную скорость на бесконечности как скорость набегающего потока. При численном интегрировании это условие ставится на некотором достаточно большом, но конечном значении координаты h, h=hinf, таком, что при дальнейшем его увеличении никаких заметных изменений профиля скорости не наблюдается. Этого легко достичь при нескольких пробных численных экспериментах (или даже так формализовать выбор hinf, чтобы вычислительная программа сама определяла нужное значение внешней границы hinf, исходя из заданной точности расчетов).

Целью дальнейших расчетов будет построение поля течения вблизи проницаемой стенки при различных значениях параметра fw - отрицательных, соответствующих вдуву (или испарению),  и положительных, как при отсосе (или конденсации).

Ограничимся далее исследованием задачи с постоянной скоростью внешнего потока  (с нулевым продольным градиентом давления) и примем соответственно нулевое значение параметра b, b=0.

В результате численного интегрирования системы ( 3) будут получены функции  F0(h), F1(h), F2(h), после чего может быть рассчитано поле течения в исходных, физических координатах, а также определено трение на стенке. Связь между физическими переменными  (т.е.  продольной  u(x, y) и поперечной  v(x, y)  составляющими скорости,  функцией тока  y(x,y)),  с одной стороны,  и автомодельными переменными Фолкнера-Скэн (т.е. безразмерной функцией тока f(h)  и  безразмерным расстоянием от стенки h) – с другой,  задается следующими уравнениями:

для продольной скорости, отнесенной к скорости внешнего потока 

( 6)

для поперечной скорости, отнесенной к скорости внешнего потока

( 7)

Уравнение ( 7), записанное для нулевого значения h, определяет связь между значениями функции тока  и поперечной скорости (скорости вдува или отсоса)  непосредственно на стенке:

 

( 8)

Вторая производная, вычисленная непосредственно на стенке, дает значение коэффициента трения, как это видно из следующего уравнения:

( 9)

Заметим, что соотношения ( 7) и ( 8) соответственно дают основание назвать величину fw – параметром проницаемости, а величину fw² - параметром трения.

Автомодельная переменная, интерпретируемая  как отношение расстояния от стенки   к толщине пограничного слоя, вычисляется через физические координаты по формуле:

( 10)

 Безразмерная функция тока f  (или F0) связана с  обычной функцией тока y уравнением:

( 11)

Физический смысл функции тока следует из определяющего соотношения

( 12)

где интеграл берется при фиксированном значении продольной координаты x, т.е. по поперечному сечению потока. Разность значений на двух линиях тока равна объемному расходу (Рис. 5).

Число Рейнольдса определяется по длине обтекаемой поверхности  L: Re = u¥L / n. Физические координаты (x, y) приводятся к безразмерному виду  отнесением к длине: X=x/L, Y=y/L.

 

Сведение краевой задачи к начальной задаче методом   sbval

Граничные условия для уравнения Фолкнера-Скэн заданы на концах области интегрирования – на стенке, при h=0,  и вдали от нее,  при ¥,  практически на некотором большом, но конечном  расстоянии h=hinf.

 Другими словами, необходимо решать краевую, или, как иногда говорят, двухточечную задачу. Численные алгоритмы для таких задач более сложны,  чем для начальных задач, когда все граничные условия заданы в начальной точке.

Применительно к начальным задачам, вычислительная математика располагает эффективным арсеналом быстрых и точных численных методов, среди которых наиболее распространены методы Рунге-Кутта с фиксированным или адаптирующимся шагом. Можно использовать этот арсенал, сводя решение краевой задачи к интегрированию серии начальных задач с пробными значениями недостающих начальных условий. Насколько удачна или неудачна каждая такая попытка, можно судить по «расстоянию», на котором оказывается вычисленное значение от предписанного граничным условием на  другом конце отрезка интегрирования.

Обсудим более конкретно эту идею применительно к нашей задаче. В начальной точке поставлены два условия – заданы значения для F0(0)  и F1(0). Недостает значения для F2(0) в начальной точке, но задано F1(hinf)=1 на конце отрезка интегрирования. Решим задачу с начальными данными, взяв в качестве недостающего условия для F2(0)  некоторое число x, может быть, просто наугад.  В результате на конце отрезка  получится решение с некоторым значением F1, которое мы обозначим как F1,inf .  Мало вероятно, что после интегрирования получится правильное единичное значение F1,inf=1, но можно считать результат  F1,inf функцией от пробного значения x . Теперь ясно, как действовать дальше. Имеется два варианта:

q              Решать уравнение F1,inf(x)-1=0 каким либо численным методом, например,  методом секущих,  понимая, что пробные значения функции F1,inf(x) находятся не просто подстановкой в какую-то формулу, но посредством обращения к процедуре численного интегрирования системы дифференциальных уравнений.

q              Стремиться минимизировать невязку (F1,inf(x)-1)2 ® 0, добиваясь нулевого значения с помощью какого-либо алгоритма оптимизации (минимизации), например, методом покоординатного спуска, или симплексным методом Нелдера-Мида, или каким-либо другим имеющимся в распоряжении методом.

В Mathcad’е имеется специальная функция sbval, посредством которой можно определить недостающие условия в начальной точке. Мы не знаем точно внутреннего устройства этой функции, но по-видимому, она работает по какому-либо из описанных выше вариантов. Как пользователям,  нам достаточно  этого общего представления, но мы должны уметь написать обращение к указанной встроенной функции (Рис. 2). Как только недостающее  начальное значение будет возвращено функцией sbval,  получившуюся начальную задачу окончательно  решают численным методом (например,  методом Рунге-Кутта с фиксированным шагом rkfixed).

Далее, воспользовавшись приведенными выше формулами  ( 6)-( 11) перехода к физическим переменным и записанной далее программой интерполяции (Рис. 4), строят векторное поле скорости и функцию тока, чтобы сделать видимым влияние вдува или отсоса на гидродинамику потока.

Набор числовых параметров для примера решения уравнения Фолкнера-Скэн будет следующим:

·               Числовой параметр, определяющий ускорение внешнего потока: b=0

·               Числовой параметр, задающий вдув  (fw<0) или отсос (fw>0):     fw=0 (это случай непроницаемой стенки)

·               Фиксированное нулевое значение координаты на обтекаемой стенке: hw=0

·               Числовой параметр, задающий внешнюю границу пограничного слоя: hinf=6

Обращение к  встроенной функции sbval показано на Рис. 2 . Последовательность действий такова. Вводятся значения параметров, записывается правая часть ( 4) системы дифференциальных уравнений ( 3). Пользовательская функция SetInit формирует вектор начальных условий в соответствии с ( 5). Пробное значение второй производной при h=0 указывается как компонент вектора с нулевым индексом x0 ( в других задачах недостающих начальных условий может быть больше, чем одно, поэтому  x должно быть вектором). Пользовательская функция discrepancy формирует невязку на конце отрезка интегрирования ((F1-1) должно быть равно нулю, см. ( 5)). Далее записывается обращение к функции sbval, а результат, т.е. недостающее начальное условие, помещается в вектор MissingInitCond. Напомним, что это недостающее условие есть значение второй производной на стенке, т.е.  F2(0).

 

 

Рис. 2. Применение функции sbval, возвращающей недостающие условия в начальной точке (F2(0)=0.4696 для непроницаемой стенки).

 

Решение начальной задачи методом rkfixed

Поскольку теперь известны все начальные условия, можно обратиться к какой-либо встроенной функции численного интегрирования, например rkfixed (Рис. 3). Параметрами этой функции являются: вектор начальных условий InitCond, координаты начальной и конечной точек, число шагов интегрирования N и вектор-функция правой части системы дифференциальных уравнений. Заметим, что фрагмент вычислений на Рис. 3 является продолжением документа Mathcad на Рис. 2, где уже была представлена правая часть D.

 

 

Рис. 3. Распределение продольной скорости в пограничном слое на непроницаемой поверхности

Результат интегрирования возвращается как массив S, составляющими которого являются векторы независимой переменной,  безразмерных функции тока, продольной составляющей скорости и напряжения трения. Для наглядности, профиль продольной скорости построен таким образом, что обтекаемую стенку следует отождествить с горизонтальной осью координат, а поток считать движущимся слева направо. По вертикальной  оси отложены значения  безразмерного расстояния от стенки h. Значения продольной скорости  отсчитываются по горизонтальной оси. В действительности построен график типа  Vector Field Plot (векторное поле), но горизонтальные стрелки продольной составляющей скорости нарисованы так часто, что сливаются, образуя эпюру скорости. Последние три строчки перед диаграммой на Рис. 3 подготавливают массивы данных для построения такой диаграммы.

Итак, на Рис. 3 представлен классический профиль скорости в пограничном слое на непроницаемой стенке.  Далее будут проведены расчеты для проницаемой поверхности, и мы увидим, как вдув или отсос управляют профилем скорости в пограничном слое. Но предварительно представим результаты еще в двух формах (Рис. 5). Во-первых, построим наглядную картину развития векторного поля течения вдоль обтекаемой поверхности, и, во-вторых, построим линии тока. 

 

Построение поля течения

Для построения поля течения в исходных, физических координатах необходимо заполнить область течения квадратной координатной сеткой (x, y), для каждого узла сетки вычислить по формуле ( 10) значения автомодельной переменной, а затем по формулам ( 6), ( 7) и ( 11) найти значения продольной и поперечной проекций вектора скорости, а также функции тока y. Эти вычисления реализованы в Mathcad-программе на  Рис. 4.

 

 

Рис. 4. Программа интерполяции для расчета поля течения в физических координатах.

 

Результаты расчетов представлены на Рис. 5. Поскольку на передней кромке обтекаемой пластины (при X=0) решение имеет особенность, построение начинается с некоторого малого конечного значения Xmin.

На верхнем рисунке видно, как постепенно по ходу потока утолщается пограничный слой, внутри которого происходит все изменение скорости – от нуля на стенке до скорости набегающего потока. Линии тока (нижний рисунок) идут почти параллельно обтекаемой поверхности, хотя заметно их некоторое отклонение вверх по ходу потока вследствие  подтормаживания,  «налипания» жидкости на твердой стенке.

И продольная, и поперечная составляющие скорости обращаются в ноль на самой твердой непроницаемой поверхности.

Напряжение трения, пропорциональное производной (u/y), достигает максимального значения на стенке и обращается в ноль на внешней границе пограничного слоя. В соответствии с уравнением  ( 9),  в котором  следует сделать замену (см. значение переменной MissingInitCond на Рис. 2)

,

коэффициент сопротивления на непроницаемой пластине выражается формулой:

 

( 13)

 

 

 

Рис. 5. Векторное поле течения и функция тока для непроницаемой стенки.

 

Пограничный слой на проницаемой поверхности.

При вдуве через проницаемую стенку (см. Рис. 6, Рис. 7, Рис. 8) наблюдается резкое изменение структуры течения в пограничном слое. Профиль скорости (Рис. 6) приобретает характерную S-образную форму. Слой сдвига, в котором происходит резкое изменение скорости и где касательные напряжения трения отличны от нуля (т.е. собственно пограничный слой), оказывается оттесненным от  обтекаемой поверхности. На самой же стенке напряжения трения уменьшается и в пределе, при так называемом критическом вдуве, обращается в ноль. Значение параметра проницаемости  fw= -0.7 , принятое в расчете, близко к критическому, а значение коэффициента трения практически нулевое (f w²=F2(0)=0.05458 , см. Рис. 6, сравнить со значением  f w²=F2(0)=0.4696 для непроницаемой стенки).

Чтобы сделать структуру течения вблизи стенки более отчетливой, на Рис. 8 показано векторное поле скорости в увеличенном масштабе. Следует отметить, что продольная составляющая на самой стенке по-прежнему равна нулю, т.е. выполняется условие прилипания. Однако  поперечная составляющая, в отличие от непроницаемой поверхности, уже ненулевая. Поэтому на самой стенке вектор скорости ненулевой и направлен по нормали к поверхности.

 

 

Рис. 6. Распределение продольной скорости в пограничном слое при вдуве

( fw=-0.7,   f w²=F2(0)=0.05458 )

 

 

 

Рис. 7. Векторное поле скорости и функция тока при вдуве.

 

 

 

Рис. 8. Векторное поле скорости вблизи стенки при вдуве.

 

Течение в пограничном слое с отсосом показано на Рис. 9, Рис. 10, Рис. 11. Мы вновь замечаем существенные отличия от непроницаемой стенки. Теперь пограничный слой прижимается к стенке, его толщина заметно меньше, чем для непроницаемой поверхности и, тем более, для течения со вдувом.

Благодаря большей крутизне профиля скорости, значительно возрастает трение на проницаемой поверхности в случае отсоса. Вторая производная на стенке, пропорциональная трению, принимает значение fw²=F2(0)=7  (при fw= 7). Это очень большая величина, если сравнить со значением  f w²=F2(0)=0.4696 для непроницаемой стенки.

На Рис. 11 показано  векторное поле скорости в увеличенном масштабе для случая течения с отсосом. Как всегда, благодаря условию прилипания продольная составляющая на самой стенке равна нулю. Поперечная составляющая, в отличие от непроницаемой поверхности, ненулевая и принимает отрицательные значения.

 

 

Рис. 9. Распределение продольной скорости в пограничном слое с отсосом

( fw=-7,   fw²=F2(0)=7.0692 ).

 

 

Рис. 10. Поле скорости и функция тока при сильном отсосе.

 

 

 

Рис. 11. Поле скорости вблизи стенки при сильном отсосе.

 

В расчетах, проведенных выше для обтекания непроницаемой поверхности и для случаев со вдувом и отсосом через стенку, варьируемым параметром было значение функции тока на стенке fw, которое мы называли параметром проницаемости (см. уравнение ( 8)). При численном интегрировании краевой задачи каждый раз для заданного параметра fw получалось соответствующее значение второй производной на стенке fw², которое является мерой трения (см. уравнение ( 9)). Результаты таких вычислений представлены графически на Рис. 12 в форме зависимости fw²(fw).

Как видно из графика, для больших положительных значений параметра проницаемости, т.е. для сильного отсоса, обнаруживается следующая асимптотика: fw² Þ fw.

 

 

Рис. 12. Зависимость параметра трения от поперечного потока

 

Из уравнений ( 8) и ( 9) следует, что в указанном асимптотическом пределе сильного отсоса

( 14)

или, что то же самое

( 15)

Физический смысл асимптотической формулы ( 15) для трения можно пояснить следующим образом. Заключенная в квадратные скобки часть этого выражения есть плотность поперечного потока массы, поступающего из внешнего потока в проницаемую стенку, а второй сомножитель есть продольный импульс, «содержащийся» в каждой единице массы этого поперечного потока. Произведение этих величин составляет плотность потока импульса через поверхность стенки, что равносильно заданию трения на стенке.

Замечательным свойством предельного закона ( 14) или ( 15) является независимость напряжения трения (или коэффициента трения) от вязкости жидкости (сравните асимптотическую  формулу ( 14) с формулой для вязкого трения на непроницаемой стенке  ( 13)).

Другой важный результат уже обсуждался выше – это обращение в ноль трения на поверхности при сильном вдуве. На графике (Рис. 12) этот случай оттеснения (или отрыва) пограничного слоя виден в области предельных отрицательных значениях параметра вдува  fw, где параметр трения fw² обращается в ноль.

 

Уравнение теплового пограничного слоя

Если температуры стенки и потока жидкости (Рис. 1) неодинаковы, то возникает тепловой поток, который в инженерных расчетах определяют по уравнению Ньютона-Рихмана

( 16)

где a - коэффициент теплоотдачи, мера интенсивности процесса конвективного теплообмена между твердой поверхностью и потоком жидкости. Разность температур стенки и набегающего потока называют температурным напором.

Коэффициент теплоотдачи является сложной функцией скорости и режима течения, формы и размеров обтекаемой поверхности, теплофизических свойств жидкости. Можно сказать, что разработка методов расчета коэффициента теплоотдачи – это главная задача теории конвективного теплообмена.

Подобно тому как уравнения ( 1), ( 2) описывают поле течения в гидродинамическом пограничном слое, аналогичное по форме уравнение  ( 17) с краевыми условиями  ( 18) определяет температурного поле в тепловом пограничном слое вблизи обтекаемой поверхности:

 ( 17)

 ( 18)

Автомодельная переменная имеет прежний смысл (см. ( 10)), искомая функция g представляет собой безразмерную температуру в потоке жидкости:

,

( 19)

а числовой параметр  Pr  (число Прандтля) – есть отношение кинематической вязкости к коэффициенту температуропроводности жидкости.

Несложные вычисления показывают, что коэффициент теплоотдачи a  и его безразмерная форма – число Нуссельта Nu -  выражаются через производную от безразмерного профиля температуры на стенке  g`w следующими формулами, которыми можно воспользоваться после интегрирования краевой задачи  ( 17),  ( 18):

 

( 20)

Как видно из  ( 17),  для интегрирования необходимо располагать решением  f(h) для поля течения. В этом разделе мы ограничимся исследованием теплообмена при обтекании непроницаемой стенки. Поскольку непосредственно необходима только функция тока f(h), применим  для интегрирования задачи  ( 1), ( 2)  встроенную функцию Odesolve, как показано на Рис. 13.

 

Рис. 13. Гидродинамический пограничный слой на непроницаемой стенке.

 

Собственно обращение к решателю Odesolve занимает всего одну строчку текста. Поскольку полученное распределение f в дальнейшем понадобится многократно при вариантных расчетах, целесообразно сохранить результаты интегрирования в форме сплайн-функции fspl . Чтобы более отчетливо представить область пограничного слоя, построено также распределение скорости и отмечено расстояние от стенки h=3.472, на котором завершается 99% всего изменения продольной скорости. Для определения этой границы попутно понадобилось решить нелинейное уравнение методом root. За пределами слоя h=6 скорость уже практически постоянна, а функция тока линейно увеличивается с расстоянием от стенки. Поэтому хорошим приближением для функции тока при любых h будет пользовательская функция fa, которая и используется далее при интегрировании уравнения теплового пограничного слоя.

Решение задачи о температурном поле представлено на Рис. 14. Ориентация осей координат на графике выбрана так, что обтекаемая поверхность совпадает с горизонтальной осью, вдоль которой откладываются значения безразмерной избыточной температуры в потоке. По вертикальной оси отложено безразмерное расстояние от стенки.

Расчеты проведены в широком диапазоне чисел Прандтля. Отметим, что близкие к единице значения числа Прандтля характерны для газовых теплоносителей, много большие единицы – для вязких малотеплопроводных органических жидкостей, много меньшие единицы – для жидких металлов. График на Рис. 14 показывает, что  разделение теплоносителей на такие три группы имеет отчетливый физический смысл.

В случае Pr = 1 профили скорости и температуры полностью совпадают  и толщины гидродинамического и теплового пограничных слоев оказываются одинаковыми.

При Pr << 1 (см. числовой пример с Pr = 0.1) толщина теплового слоя намного больше толщины гидродинамического. Можно сказать, что течение в основных пределах теплового пограничного слоя происходит с максимальной равномерной скоростью, за исключением самой близкой окрестности стенки.

При Pr >> 1 (см. числовой пример с Pr = 10) толщина гидродинамического слоя намного больше толщины теплового. Тепловой слой находится на самом дне гидродинамического слоя. Можно сказать, что течение в пределах теплового пограничного слоя происходит с относительно малой скоростью, уменьшенной из-за близости к стенке, непосредственно на которой скорость обращается в ноль.

 

Рис. 14. Температурное поле при различных значениях числа Прандтля.

 

Рис. 15. Безразмерный  градиент температуры как функция от числа Прандтля

Интегрирование уравнения теплового пограничного слоя проводилось многократно для различных значений числа Прандтля. В качестве основного результата была получена таблица на Рис. 15, в верхней строке которой указаны числа Прандтля, а в нижней – соответствующие значения безразмерного градиента температуры g`w.

Эти результаты представлены также графически, причем сопоставлены данные численного интегрирования, асимптотические выражения для малых и больших значений числа Прандтля, а также приближенное универсальное уравнение g`appr для градиента на стенке g`w, которое получено методом интерполяции между асимптотами и пригодно во всем диапазоне чисел Прандтля. Значения g`w подставляются в формулы ( 20), чтобы вычислить искомое значение коэффициента теплоотдачи, а затем и плотности теплового потока по уравнению ( 16).

Полезно по результатам  теории получить представление о реальных цифровых значениях характерных величин для трех теплоносителей, представляющих три группы, указанные выше, а именно, для ртути, воздуха и масла МК, а также для самого важного теплоносителя – воды  (

Рис. 16).  Теплофизические свойства взяты для нормальных условий. Длина обтекаемой стенки и скорость были выбраны небольшими, чтобы выполнялось условие ламинарного режима. Коэффициент теплоотдачи рассчитывался по формуле ( 20):

.

 

Рис. 16. Теплоотдача к различным теплоносителям.

 

Результаты вычислений показывают, как сильно зависит интенсивность теплообмена от физических свойств жидкости. Теплоотдача к потоку ртути (жидкого металла с высокой теплопроводностью) на три порядка превышает теплоотдачу к воздуху. Видно, что и вода является гораздо лучшим охладителем, чем вязкое масло или воздух с его малой плотностью и теплопроводностью.

Другой важный вывод теории – существенный рост интенсивности теплообмена при увеличении скорости течения, пропорционально корню квадратному из значения скорости потока (см. формулу ( 20)).

Хотя эти результаты получены для ламинарного режима, качественный характер полученных зависимостей сохраняется и для более важных для практики ( и более сложных для расчета) турбулентных потоков.

 

 

Неприятности при решении краевой задачи посредством функции Odesolve

Функция Odesolve, появившаяся в версии Mathcad 2000, предназначена для интегрирования обыкновенного дифференциального уравнения (не системы уравнений), в том числе с двухточечными граничными условиями. Обращение к этой функции по форме записи практически не отличается от обычной математической нотации и выглядит существенно проще, чем при использовании способа sbval. Однако к применению этой функции следует относиться с осторожностью, как показывает пример на Рис. 17.

Решение уравнения Фолкнера-Скэн записывается в одну строчку. При обращении к Odesolve  действительно возвращается функция, а не массив, как в случае других решателей (rkfixed и др.). Проблема состоит в том, при обращении к этой функции могут получиться неправильные результаты.

Первая ошибка выявляется при вычислении первой и второй производных в начальной точке (строка 2 на Рис. 17). Видно, что не воспроизводится нулевое значение первой производной, заданное начальным условием. Для второй производной получается ноль вместо правильного значения 0.4696 (см. Рис. 3). Это выглядит уже не как погрешность вычислений, а как грубая ошибка.

Третья строчка показывает, что возвращаемый объект f не является обычной функций, т.к. не удается применить стандартную операцию разложения в ряд.

Фрагменты 4-6 иллюстрируют странности, выявляющиеся при попытках последовательно применить сплайновую интерполяцию и уже затем дифференцирование. Фрагмент 4 показывает, что при малом шаге интерполяции получаются такие же неправильные результаты, как при прямом вычислении производных в строке 2.

Фрагмент 5 показывает, что можно подобрать шаг интерполяции, при котором результаты будут удовлетворительными. Однако следующий фрагмент 6 опять вызывает недоумение: незначительное изменение шага приводит к резкому изменению результата.

По-видимому, для функции Odesolve нужна «заплатка» от производителей.

 

 

Рис. 17. Функция Odesolve

 

Заключение

Картины векторного поля течения наглядно показывают, что при заданном в одном из примеров большом отрицательном значении параметра вдува  fw = -0.7  имеет место явление оттеснения пограничного слоя. В непосредственной близости от стенки вектор скорости направлен практически вертикально, и только на некотором удалении от нее развивается продольное течение со сдвигом. Такие режимы со вдувом теплоносителя через стенку используются для защиты поверхностей от воздействия высокотемпературного потока газа или для защиты от химически агрессивных сред.

Варьируя в расчетах значение параметра fw, можно наблюдать  влияние поперечного потока массы на гидродинамику пристенного течения. При fw=0 получается стандартная картина обтекания непроницаемой стенки, а при положительных значениях - картина течений с отсосом . Заметим, что в аэродинамике отсос используется для предотвращения опасного явления отрыва пограничного слоя при обтекании крыла с большими углами атаки, как при посадке самолета.

Интересные эффекты можно  наблюдать при совместном влиянии параметров ускорения b и проницаемости fw. Варьируя эти два параметра, инженер-проектировщик может осуществлять целевое управление пограничным слоем, например, применяя отсос, чтобы предотвратить отрыв при течении против градиента давления.

Решения уравнения  Фолкнера-Скэн и аналогичного по структуре уравнения для теплового пограничного слоя выявляют основные закономерности трения и теплообмена, показывая, какую роль в этих процессах играют физические свойства жидкости, размеры обтекаемой поверхности, скорость течения. Специальным образом представленные в форме  законов трения, теплообмена и массообмена для ламинарных течений,  результаты решения используются в рамках приближенных интегральных методов, ориентированных на инженерные приложения [4, 5] .

 

Литература

1.      Falkner V.M., Skan S.W.. Some approximate solutions of the boundary layer equations. Phil. Mag. 12, 865-896 (1931)

2.      Spalding D.B., Evans H.L.. Mass Transfer throgh Laminar Boundary Layers. Int.J. Heat Mass Transfer. Vol.2.  pp. 199-221. Pergamon Press 1961.

3.      Na T.Y.. Computational Methods in Engineering Boyndary Value Problems. Academic Press. 1979. (русский перевод: На.Ц.. Вычислительные методы решения прикладных граничных задач. Москва. Мир. 1982. 296 с.)

4.      Кутателадзе С.С., Леонтьев А.И. Тепломассообмен и трение в турбулентном пограничном слое. М.: Энергоатомиздат, 1985. 320 с.

5.      Солодов А.П. Интегральный метод решения задач пограничного слоя. Москва. Изд-во МЭИ. 1992. 79 с.